The Intrinsic Shape Distribution of a Sample of ElUptical Galaxies 



Jakob Bak^ and Thomas S. Statler^ 
Department of Physics and Astronomy, Ohio University, Athens, OH 45701 

o ■ 

Q ■ Received ; accepted 

(N 



> 
oo 

m 
o 
o 
o 

o 



^bak@hehos. phy.ohiou.edu 
^tss@hehos. phy.ohiou.edu 



- 2 - 



ABSTRACT 



We apply the dynamical modeling approach of Statler (1994b) to 13 elliptical 
galaxies from the Davies and Birkinshaw (1988) sample of radio galaxies to 
derive constraints on their intrinsic shapes and orientations. We develop an 
iterative Bayesian algorithm to combine these results to estimate the parent 
shape distribution from which the sample was drawn, under the assumption 
that this parent distribution has no preferred orientation. In the process we 
obtain improved estimates for the shapes of individual objects. The parent 
shape distribution shows a tendency toward bimodality, with peaks at the oblate 
and prolate limits. Under minimal assumptions about the galaxies' internal 
dynamics, 35% of the objects would be strongly triaxial (0.2 < T < 0.8). 
However, the parent distribution is sensitive to the assumed orbit populations in 
the galaxies. Dynamical configurations in which all galaxies rotate purely about 
either their long or short axes can be ruled out because they would require the 
sample to have a strong orientation bias. Configurations in which the mean 
motion about the short or long axis is either "disklike" — dropping off away from 
the symmetry planes — or "spheroidlike" — remaining roughly constant at a given 
radius — are equally viable. Spheroidlike rotation in the long-axis or short-axis 
tube orbits significantly lowers the abundance of prolate or oblate galaxies, 
respectively. If rotation in ellipticals is generally disklike, then triaxiality is rare; 
if spheroidlike, triaxiality is common. 



Subject headings: galaxies: elliptical and lenticular, cD — galaxies: kinematics 
and dynamics — galaxies: structure 



-3- 



1. Introduction 

Over the years a number of attempts have been made to derive the intrinsic shape 
distribution of eUiptical galaxies from observations (Hubble 1926, Sandage et al. 1970, 
Noerdhnger 1979, Marchant & Olson 1979, Richstone 1979, Binggeh 1980, Binney & de 
Vaucouleurs 1981, Olson & de Vaucouleurs 1981; for a review see Statler 1996). Generally 
the results of these efforts have been ambiguous, and interest in the problem waned 
somewhat in the 1980s. But more recent developments have sparked renewed attempts to 
crack this classic chestnut. Among these developments are the recognition that halo shapes 
may serve as a diagnostic of galaxy formation physics (Dubinski & Carlberg 1991, Weil & 
Hernquist 1996), and indications that Hamiltonian chaos, dissipation, or both may either 
force triaxial equilibrium configurations to evolve slowly toward axisymmetry or render 
them altogether impossible (Dubinski 1994, Merritt & Pridman 1996, Merritt & Quinlan 
1998). 

Studies of central surface brightness profiles using HST suggest that the fundamental 
properties of elliptical galaxies may be bimodally distributed. There appears to be a 
dichotomy between high- luminosity, slowly rotating systems with shallow central cusps and 
boxy isophotes, and lower luminosity, rotationally supported systems with steeper cusps 
and a tendency for diskiness (Lauer et al. 1995; see also Kormendy & Bender 1996). Since 
rapid rotation and triaxiality are generally regarded as being incompatible, one might 
anticipate a bimodal distribution of triaxialities. Tremblay & Merritt (1996) find that low- 
and high-luminosity elliptical galaxies have different distributions of apparent ellipticity, 
which would imply different distributions of true shapes. Merritt & Tremblay's work joins 
that of Fasano & Vio (1991), Ryden (1992, 1996), and Fasano (1996) as successors to the 
classical photometric approaches pioneered by Hubble, Sandage, and others. 

However, photometric methods, while effective in constraining the distribution of 
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overall flattenings, reveal little about the frequency of axisymmetry vs. triaxiality in the 
population. Uncovering this information requires the use of kinematic data and dynamical 
models to connect the kinematics to the shape of the gravitational potential. In the rare 
cases where well-defined, equilibrium gas disks are present, emission-line kinematics can 
yield excellent constraints on the shape of the potential if one assumes that the gas is on 
closed orbits (Bertola et al. 1991). But for the majority of ellipticals, methods relying 
primarily on stellar kinematics are essential. Approaches of this type were originated by 
Binney (1985) and enlarged upon by Pranx et al. (1991) and Tenjes et al. (1993). Statler 
(1994a, 1994b) introduced major refinements, including improved dynamical models and a 
Bayesian approach to model fitting. This method has great potential to place quite narrow 
constraints on the triaxialities of individual galaxies for which very high quality stellar 
kinematic data are available. Unfortunately, the number of such galaxies is still very small, 
and is likely to increase at only a modest pace in the short term. 

Our goal in this paper is to see what can be learned from the larger sample of galaxies 
with stellar kinematic data of less-than-ideal quality already in the literature. We focus 
on the Davies & Birkinshaw (1988, hereafter DB) sample of radio ellipticals, all of which 
have kinematic data on multiple position angles and are photometrically well studied. In 
the process, we extend the statistical methods of Statler (1994b) and show how to estimate 
the parent shape distribution from a sample of galaxies for which the data may be very 
inhomogeneous. Our method will thus continue to be generally applicable as new data are 
obtained. 

In the next section of the paper we describe the general statistical approach for 
determining the parent distribution of a set of intrinsic quantities from measurements of 

related, but different, observable quantities, and show the particular application of this 
approach to the shape problem. In § 3, we discuss our treatment of the data and define a 
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subsample of the DB galaxies which we are able to model reliably. Section 4 presents the 
results for the parent distribution, and examines systematic effects relating to unknown 
aspects of the stellar dynamics. Section 5 compares our results to those of previous studies, 
and § 6 sums up. 



2. Estimating the Parent Distribution 
2.1. Basic Idea 

In previous papers (Statler 1994b, 1994c, Statler et al. 1999) we describe a Bayesian 
approach to inferring the intrinsic shapes of individual elliptical galaxies, based on 
dynamical models that predict the mean radial velocity field (VF) by solving the equation 
of continuity for the stellar "fluid." For triaxial systems with negligible figure rotation, the 
streamlines of the mean motion in the main families of circulating orbits are dictated by the 
triaxiality of the mass distribution; thus for a given shape, orientation, and set of boundary 
condition parameters describing the internal orbit populations, the line-of-sight VF can be 
calculated. The calculation of the VF is very fast, so the multidimensional parameter space 
can be adequately explored. 

For each set of parameters, we calculate the probability that the observed VF and 
surface brightness distribution would, with the known observational errors, be obtained from 
the corresponding model. This yields a multidimensional likelihood function L(T, c, f2,d), 
where T is the triaxiality of the total mass distribution, c is the short-to-long axis ratio of 
the luminosity distribution, Q = {6, (p) is the orientation of the galaxy, and the vector d 



represents the remaining dynamical parameters (see § |2.4|) . The likelihood is integrated 



over an assumed prior distribution in the dynamical parameters d and an isotropic prior 
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distribution in Q to give a two-dimensional likelihood L(T, c).0 To obtain the Bayesian 
estimate of the galaxy's shape, L{T, c) is multiplied by a model for the parent shape 
distribution and normalized. In most of our previous work, this model has been a flat 
distribution, F{T, c) = const, meaning that we have estimated the shape of each galaxy 
in isolation. The likelihood Lj(T, c) for each galaxy i is basically a data point with an 
error ellipse, i.e., a measurement of T and c. Of course, in this case the errors are strongly 
non-Gaussian since we work with the actual probability distribution. The goal now is 
to combine these measurements for a sample of galaxies into an estimate of the parent 
distribution from which the sample was drawn. 

Our algorithm is conceptually simple. We start with a fiat model for the parent 
distribution, F(T,c) = const. The Bayesian posterior probability density Pi(T,c) for each 
galaxy is the normalized product of Lj with the model parent F. (At this stage, = L,.) 
We stack the Pj's on top of each other, add them up, smooth the sum with a nonparametric 
smoothing spline, and normalize the result. This gives us an improved model for the 
parent distribution F, which we multiply by the Lj's and feed iteratively into the same 
procedure. Note that, after the first iteration, the statistical estimate of the shape of each 
galaxy in the context of the whole sample. Pi, is different from the estimate of its shape in 
isolation, Lj. This difference arises from the requirement that the sample be drawn from an 
isotropic distribution of orientations. Note also that since all operations are performed on 
distributions that are already integrated over fl, the isotropy of the parent distribution is 
guaranteed. 

A one dimensional toy problem can demonstrate that this algorithm works, even when 
the Li distributions are strongly non-Gaussian. Consider a set of objects, each with a value 
of some intrinsic property X between and 1. X is not measurable, but a related quantity 



■^See Statler et al. (1999) for a more rigorous formulation of this statement. 
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X is. Suppose that for any object an observer has a uniform probabihty of measuring any 
value of X between and X. It is easy to show that a single measurement Xi implies a 
likelihood Li{X) that is zero for X < Xi and proportional to X~^ for X > Xi. Figure |I] 
shows the algorithm in action. The likelihoods Li{X) for 9 measured x values are shown in 
Fig. |l]a. These functions are multiplied by the initially fiat parent distribution (b), summed 
(c), and smoothed to produce the new parent (d). The result after 20 iterations (e, solid 
line) is a decent representation of the true parent distribution (dotted line). The functions 
Pi{X) giving the estimates of the X values for the individual objects (f) differ from the 
original Lj's but are consistent with the parent distribution. 

2.2. Statistical Rationale 

The stack-smooth-iterate algorithm is a general technique that is closely related 
to Lucy's Method (Lucy 1974) and penalized likelihood (Wahba & Wendelberger 1980, 
Silverman 1986, Green & Silverman 1994) Readers uninterested in the statistical details 
are welcome to skip directly to § 2.5. 

Imagine a population of objects, each of which possesses some value of an intrinsic 
property (or set of properties) X, distributed according to the parent distribution Fp{X). 
Let there be an observable quantity (or set of quantities) x which is related to X by a 
conditional probability distribution P{x\X). The distributions are normalized so that 
/ dX Fp{X) = J dx P{x\X) = 1. /^From models, we calculate the likelihood L{X\x) that 
a given measurement of x was obtained from an intrinsic value X. We assume that the 

^Our method is similar, but not identical, to a well-established approach developed by 
Wahba & Wendelberger (1980). We are grateful to the referee for directing us to this 
important paper. 
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likelihoods are explicitly normalized so that / dX L{X\x) = 1. If the models take into 
account all physical effects and measurement errors exactly, then L{X\x) and P{x\X) are 
mathematically the same function. For measurements x, the estimate of X is given by the 
posterior density, 

p,^^ _ F{X)L{X\x) 
^ ' JdXF{X)L{X\x)' ^ ' 

where F{X) is the current estimate of the parent distribution. This is the standard 
Bayesian approach for individual objects. 

The goal is to find the parent distribution F{X) that maximizes the joint probability 
of obtaining the measurements Xi for the set of n objects i = 1, . . . ,n. The logarithm of 
this probability is given by 

\nV = InJ] / dX F{X)P{xi\X) = Y.^^ j dx F{X)L{X\xi). (2) 

i i 

If we write the measurements in terms of a distribution of observables, W{x) = J2i — Xi), 
then we can interpret the integral 

Wmix) = J dXF{X)L{X\x) (3) 

as giving the distribution of observables that would be predicted if the model parent 
distribution F{X) were correct. With these definitions, equation becomes 

\nV = J dxW{x)\nWm{x). (4) 

If we subtract the constant / dxW{x) lnW^(a;) from the quantity in equation (^), we get the 
Lucy iJ-function, 

H, = JdxW{x)ln^. (5) 

In the statistics literature, —H^ is known as the "Kullback-Leibler information distance" 
between model and data (Silverman 1986). 
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Lucy's method works to increase (decrease the information distance) by iteratively 
applying the rule 

Fnew(X) = F{X) J dx ^^^L{X\x). (6) 
Using equations (|l]) and (|), this can be written as 

FneAX) = J dxW{x)P{X) = ^P,(X). (7) 

i 

Thus one iteration of Lucy's method is identical to our scheme of "stacking" the posterior 
densities. In the absence of smoothing, our approach will seek a parent distribution that 
maximizes the likelihood of the observed sample. 

For a finite sample, however, the maximum-likelihood parent distribution will be a 
set of spikes at the maximum of each L{Xi\xi), so it is ill-advised to iterate without a 
penalty function that enforces smoothness. Moreover, for a realistically small sample, 
there is not a unique F{X) that maximizes InP unless such a penalty function is present 
to lift the degeneracy. To implement this penalty at each iteration, we regard Fnew(-^), 
computed according to equation ([^), as a noisy realization of an underlying smooth function 
This function is estimated using a smoothing spline, and Fnew(-^) is replaced by 
its smoothed counterpart. 



2.3. Smoothing Splines and Cross- Validation 

Smoothing splines may be defined in any number of dimensions; here we summarize 
the two-dimensional case. This discussion is adapted from Green & Silverman (1994) and 
Silverman (1986). 

The estimate of the parent distribution, F^cw{X), is defined on a discrete grid of n 
values of X. (We drop the subscript "new" in what follows for brevity.) In our case, 
X = {T,c), and we have F{Ti,ci),..., F{Tn,Cn), from which we want to determine the 
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underlying smooth function F^{T, c). For a trial function g{T, c), a penalized sum of squares 
of the residuals is given by 

n 

S{9) = E{^(^- c,) - giT,, c,)Y + aJ{g), (8) 
1=1 

where a > is a "rate of exchange" between the usual goodness of fit measure and the 
penalizing function J{g), given by 




The penalizing function measures the rapid variation and departure from local linearity in 
g. The functions which minimize S{g) are known as thin plate splines, which are analogous 
to natural cubic splines in one dimension. Algorithms for calculating thin plate splines are 
implemented in the routines DTPSS and DPRED in the GCVPACK package (Bates et al. 
1987), which is available from Netlib.Q 

The problem of finding F^{X) is now reduced to determining a single parameter a. 
Likelihood cross validation provides a method for doing this automatically. The premise is 
that the best estimate of a should produce the distribution which best predicts all future 
data points. Since one is generally not gifted with prescience, one proceeds by removing one 
measurement from the sample and calculating the likelihood that that measurement would 
be obtained in the parent distribution found from the other measurements. Repeating this 
procedure for each measurement in turn and then averaging the likelihoods yields the cross 
validation (CV) score. 

In our case, the measurements are the individual normalized likelihoods L{Xi\xi). We 
remove the ith measurement from our data set and create a new distribution using 
the methods above. For a given value of a, the likelihood that a single L{Xi\xi) is drawn 



^http:/ /netlib.org/gcv 
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from the smoothed model parent distribution Ff_-{X; a) is given by 



In 



(y dXFli{X;a) L(X|x,)) 



(10) 



Averaging the >C_i's gives the hkehhood cross vahdation score, 



n 



CV{a) = n 



-1 



(11) 



and maximizing CV{a) provides the best estimate for a. To ensure uniqueness of the final 
result, we compute the maximum of CV{a) only once, on the first iteration, and fix the 
smoothing parameter for all subsequent iterations to its initial value. In some cases there 
is not a unique maximum in CV{a); instead, CV{a) is nearly flat up to some ckq, beyond 
which it turns over. We set a to its turnover value, and we see no indication that this 
choice biases the results. 



Our numerical implementation follows from that described in Statler (1994b). The 
treatment of individual galaxies is essentially the same, except for details noted in § 3 
below. The grid of dynamical parameters used in the models is also the same as in the 
earlier work; to aid the reader in § 4 we give a brief overview here. 

The models assume that (1) rotation of the figure (i.e., tumbling) is neghgible; (2) 
short-axis tube and long-axis tube mean motions can be represented by confocal streamlines 
(Anderson & Statler 1998); (3) the luminosity density pl is stratified on similar ellipsoids, 
PL{r,9,(f)) — pL{r)p\{9,(t)); and (4) the velocity field obeys a "similar fiow" ansatz outside 
the tangent point for a given line of sight, v(r, 9, 0) = v{r)Y*{9, (f)). The last two assumptions 
are needed for projecting the models. The results are insensitive to the accuracy of these 
assumptions as long as phir) and pL{r)v{r) decrease faster than r~^. This requirement 



2.4. Implementation 
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limits the validity of the models to regions where the rotation curve is not steeply rising. 
As a further simplification we adopt power laws for the luminosity density and the velocity 
scaling law: p ~ r"*^ and v{r) ~ r~'. The index k is determined from surface photometry, 
and we nominally adopt I — (0, ± 72), omitting the I — —^ji case when /c < 2.5. It turns 
out that the results are not very sensitive to either of these parameters. 

Remaining properties of the phase space distribution function are described by a scalar 

constant C and a function of one variable v*{t\ These parameters describe the mean 

velocity across the xz plane on one fiducial shell, which in turn determines the velocity field 

over the whole shell once the triaxiality T and the luminosity density are specified. The 

"contrast" C is defined as the ratio of the y component of the mean velocity on the x axis 

to that on the z axis, on the fiducial shell. The function gives the angular dependence 

of the mean velocity across the xz plane on the fiducial shell. The variable t is a rescaled 

polar angle, given, for spherical shells, by 

f 2-2iHii ^<sin-^VT, 
t=\ ^ ' (12) 

\ Q > sin-^ VT, 

where 9 is the usual polar angle. The relation for ellipsoidal shells is given in § Section 3.1 
of Statler (1994b). By definition, v*{0) = C and v*{2) = 1. 

The model grid comprises 8 different assumptions for the variation of C with intrinsic 
shape. In four of these C is constant: C — (long-axis tube dominated), 0.5, 1, and 
infinite (short-axis tube dominated). Four more functional forms for C(T, c) are introduced 
to mimic certain self-consistent models, and are given in equations (11) - (14) of Statler 
(1994b). The function v*{t) is taken to be either piecewise-constant or piecewise-linear 
in each of the intervals [0, 1) and (1,2] (in the linear cases dropping to zero at t — 1). 
This function describes how the mean rotation speed in each of the tube orbit families 
declines away from the symmetry plane that contains its parent orbits. For example, the 
mean rotation in the Galaxy drops with height above the disk plane as one moves into 



-13- 



the more pressure-supported halo. At the other extreme, a maximally rotating isothermal 
sphere has constant rotation speed at all latitudes. Accordingly, we refer to linear v*{t) as 
"disklike" rotation, and constant v*{t) as "spheroidlike" rotation. A model can be disklike 
or spheroidlike in either short-axis or long-axis tubes. One should avoid the impression that 
disklike rotation necessarily implies a two-component structure; in an oblate disklike model, 
the mean rotation speed 45° up from the equatorial plane is half of the in-plane value, a 
much gentler transition than in a genuine disk-halo system. 

The likelihoods Lj(T, c) for each galaxy arc computed on a 20 x 20 rectangular grid on 
the intervals < T < 1 and 0.4 < c < 1. Smoothing a function over a finite domain creates 
problems near the edges unless suitable boundary conditions are imposed to minimize this 
effect. The thin plate spline does not impose any strict boundary conditions but rather 
sees the area outside of the boundaries as lacking information. The penalizing function 
J{Li{T,c)) is therefore the only part to contribute to the penalized sum of squares of 
the residuals, forcing the function Li(T, c) to be flat outside the boundaries (see section 
2.3). In practice this has the effect of biasing Li{T, c) towards closed contours and reduced 
variability near the edges (Green & Silverman 1994). The effect however is limited and our 
results show httle if no evidence of it. 

We find that, in practice, convergence of the parent distribution can be rather slow, as 
peaks grow at the expense of valleys that sink toward zero. With our sample of only 13 
objects, iterating until a stringent convergence criterion is satisfied may be dangerous. In 
order to be conservative in our conclusions regarding the frequency of triaxiality, we stop 
iterating when the maximum fractional change in F{T, c) per iteration falls below 10%. 
Typically this occurs after about 7 iterations. 



3. Data 
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3.1. Kinematics 

All of the galaxies modeled are taken from the sample of radio ellipticals for which DB 
obtained multiple position angle rotation curve measurements. The sample contains more 
E3-E4 and fewer EO galaxies than the general population and, as DB point out, includes 
an overabundance of "unusual" objects. Where appropriate we have supplemented the DB 
data with data from Franx et al. (1989), Binney et al. (1990), Bender et al. (1994) and 
Fried & Illingworth (1994). 

The published rotation curves are first oriented to match our convention that radii west 
of north are positive. Since the models assume that the rotation curves are antisymmetric, 
we fold the profiles about the center of the galaxy to reduce the formal errors in the average 
rotation velocity. For each individual galaxy we approximate by eye the radius at which 
the rotation curve flattens and use the data outside of this in the models. At large radii the 
kinematic data become unreliable for reasons which vary from galaxy to galaxy (see section 
3.3). We therefore set an outer radius beyond which we discard the data. We average the 
data points which are left between the inner and outer radii on each PA, weighted by the 
inverse square of the published errors. The uncertainty associated with the average is taken 
to be the (l/cr^)-weighted standard deviation following Statler (1994c). The inner and 
outer radii and the adopted mean velocities are given in columns 4, 5 and 9, respectively, in 
Table |1|. 

3.2. Photometry 

With the exception of the data for NGC 4839, all of the photometry is drawn from 
Peletier et al. (1990), who tabulate the ellipticity, major axis PA and surface brightness as 
functions of radius. Similiar photometry for NGC 4839 is drawn from Joergensen et al. 
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(1992). For each galaxy the adopted major axis position angle is the average between the 
inner and outer limiting radii. The ellipticities are determined by taking the unweighted 
mean in the same interval, with the standard deviation serving as the uncertainty. The 
adopted mean ellipticities and major axis PAs are given in columns 2 and 3 of Table |l|. 

The slope of the surface brightness profile is calculated by differentiating numerically. 
The surface brightness slope is then deprojected into a volume brightness slope (k) by 
adding 1. Although this is strictly valid only for pure power-law profiles, it is fine for our 
level of approximation. For most galaxies the logarithmic slope of the surface brightness 
profile is not constant in the relevant intervals, and so two values are used that span the 
ranges of k. We compute all of the models using both values. The spanning values of k for 
all of the galaxies are in columns 6 and 7 of Table |1|. 

3.3. Notes on Individual Galaxies 

Of the 14 galaxies in the DB sample, four, NGC 1600, NGC 4374, NGC 4636 and NGC 
4839, do not show any significant rotation at DB's level of accuracy. We therefore model 
them using only their photometric data. A fifth object, NGC 4278, does show significant 
rotation but is not used. It shows a 20° isophotal twist between 20" and 60" and a drop in 
rotation velocity to zero outside of 20" that conflicts with our assumption that the rotation 
curve is flat at large radii. It could be modeled but a more sophisticated method involving 
fitting at multiple radii would be required. 

Details of how we have handled the data for the remaining sample of 13 galaxies are as 
follows: 

NGC 1600, NGC 4374, NGC 4636 and NGC 4839. Because of the lack of 
any significant rotation in these galaxies, photometric data alone is used to estimate their 
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shape likelihood distributions. The triaxialities of these four galaxies are therefore poorly 
constrained. They are included on the grounds that omitting them could bias our results 
away from strongly triaxial systems, most of which are probably slowly rotating. In each 
case, the data are averaged from the center out to the largest radius for which kinematic 
data is available (see Table |l]). The ellipticity of NGC 1600, NGC 4374 and NGC 4636 
varies by about 0.10 in this interval, but the ellipticity of NGC 4839 rises steadily from 0.20 
near the center to 0.50 at 32". 

NGC 315. The turnover radius of the rotation curve is easily located at 5" by 
inspection. At 27" on PA 40 the rotation velocity is more than 3o" from the mean. Since we 
do not know if this is a real effect, we eliminate the total of six data points outside of 20". 

NGC 741. The relatively low surface brightness of this galaxy results in large 
uncertainties in the kinematic data. Nonetheless, its six position angle measurements of the 
rotation curve make it very attractive to model. Outside of 15" some rotation appears on 
PA 10 and PA 40. Although it is not at all clear that the turnover radius of the rotation 
curve has been found, all of the data from 15" to the outermost data point at 30" are used. 

NGC 1052. The data for this galaxy is the best for any in the sample. DB present 
kinematic data on four PAs, Binney et al. (1990) on the major and minor axis and Fried 
& Illingworth (1994) on the major axis. The turnover radius of the rotation curve is easily 
seen to be at 15" for this galaxy. A velocity difference of almost 50 km/sec between the 
southern and northern parts of the galaxy outside of 37" on PA 117 and PA 164 in both 
the DB and Binney et al. (1990) data sets imply that NGC 1052 may not be antisymmetric 
outside of this radius as is assumed in our models. The total of 11 questionable points 
outside of 35" are therefore eliminated. Fried & Illingworth (1994) measure the rotation 
curve to be flat out to 40" on PA 117. 

NGC 3379. The turnover radius of the rotation curve appears at about 15" so only 
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data from this radius to the outermost data point at 34" is used. 

NGC 3665. The rotation curve is fiat from inside of 5" to close to 30", but there is 
a discontinuity at 10" where the eUipticity suddenly drops from 0.35 to almost zero. The 
ellipticity then slowly rises to approximately 0.2 at 15" outside of which it is constant. 
There is also a 10° isophotal twist in the same range. Only data outside of 15" is used. 

NGC 4261. DB provide kinematic data on four position angles out to 55". The major 
and minor axis data are supplemented with data from Bender et al. (1994), who place their 
slits 4° from those of DB. The average of the sht positions of the two papers is therefore 
used when combining the two datasets. Although this does introduce some error into the 
data it is very minor compared to the uncertainty in the velocity measurements. This 
galaxy is clearly a minor axis rotator with a turnover radius of the rotation curve at 20". 

NGC 4472. The turnover radius of the rotation curve on all the PAs is at 
approximately 25". Between 3" and 30" the ellipticity of NGC 4472 increases from 0.06 to 
0.17 but is constant outside this range, so only data points beyond 30" are used. Outside of 
60" two data points on PA 160 are more than Aa from the average so we discard the points 
outside of this radius on all position angles. 

NGC 4486. This galaxy's slow rotation makes the errors relatively large in the 
velocity measurements, but outside of 20" the data is statistically consistent with a 
flattening of the rotation curve out to the last datapoint at 60". The data between these 
radii are therefore used. 

NGC 7626. This is another slow rotator. The only significant rotation is on the minor 
axis. Outside of 20" the rotation curve seems to reverse, but the errors are so large that 
the reversal is not statistically significant. We model the data only inside 20". A turnover 
radius in the rotation curve on the major axis at approximately 5" sets the inner radius. 



- 18 - 



4. Results 

4.1. The "Maximal Ignorance" Shape Distribution 

As in previous papers, we take the result from an unweighted combination of all models 
to represent the case of "maximal ignorance," i.e., minimal assumptions as to the character 
of the internal dynamics. The parent shape distribution after seven iterations is shown 
in Figure |^. The distribution is plotted in terms of T and c such that oblate spheroids, 
prolate spheroids, and spheres lie, respectively, along the right, left, and top margins. This 
distribution is bimodal, dominated by one group of nearly oblate, moderately flattened 
systems and a second group of rounder, nearly prolate systems. The valley between 
the two peaks represents a dearth of very triaxial galaxies. The bimodality is almost 
entirely a consequence of the kinematic data; to illustrate, we show in Figure ^ the result 
obtained from photometry alone, ignoring the kinematics. As discussed in the Introduction, 
photometry is effective in constraining the overall flattening distribution but reveals little 
about triaxiality. 

A more succinct description of the frequency of triaxiality in this distribution comes 
from the one- dimensional distribution F{T), obtained by integrating Figure 0a over c. The 
result is shown in Figure |^. We somewhat arbitrarily set boundaries at T = 0.2 and 
T = 0.8 to delineate "nearly oblate," "triaxial," and "nearly prolate" regions. By this 
definition, the maximal ignorance distribution is 47% nearly oblate, 18% nearly prolate, and 
35% triaxial. If we continue iterating beyond our nominal stopping criterion, the triaxial 
fraction decreases further, so we can take this as a conservative estimate of the rarity of 
triaxial systems implied by our subsample of the DB galaxies. The result is influenced 
somewhat by the four galaxies without kinematic data; if these objects are omitted, the 
fractions change to 55% oblate, 25% prolate, and 21% triaxial. 
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We can obtain some measure of whether the DB subsample is representative of the 
eUiptical galaxy population at large by calculating the expected ellipticity distribution for a 
randomly-oriented population drawn from the inferred parent. We plot this as the smooth 
curve in Figure compared with the observed ellipticity distribution from Ryden (1992). 
The two distributions are similar, though our predicted distribution contains a slight excess 
of very round galaxies. A Kolmogorov-Smirnov (KS) test implies a 14% probability that 
the observed sample was drawn from our distribution. However, the KS probability can be 
affected by details of how the mean ellipticities are defined. The ellipticities tabulated by 
Ryden are weighted by luminosity, whereas we exclude data from the brightest parts of the 
galaxies. Applying a systematic shift as small as Ae = 0.019 to our expected distribution 
would increase the KS probability to 99%. We conclude that our maximal ignorance parent 
distribution is consistent with the ellipticities of the general population of elliptical galaxies. 

The final posterior densities describing the shapes of the individual galaxies in the 
sample with rotation data are shown in Figure ^. Some well-known objects are found to 
have well constrained triaxialities; NGC 1052, NGC 3379, and NGC 4472 are probably 
oblate or nearly so. The famous minor-axis rotator NGC 4261, not surprisingly, turns out 
to be most likely prolate, though there are oblate models not excluded at the 2a level. The 
shapes of other objects are not as well constrained, and bimodal posterior densities are less 
a consequence of the kinematic data for the individual galaxies than a reflection of the 
parent distribution. 

4.2. Dynamical Configurations That Can Be Ruled Out 

Just as the marginal posterior densities describing the shape of each galaxy (Fig. 
H) can be computed for a given parent distribution, we can compute marginal densities 
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describing the orientation of each galaxy according to 

Pi{n) = j dT j dc ^F{T, c)U{T, c, (13) 

The 4-dimensional hkehhoods Li{T, c, VL) are obtained by integrating the original likelihood 
function Lj(T, c, fi, d) over the dynamical parameters. The factor l/47r reflects the assumed 
isotropy of the parent distribution; in other words we have assumed that the 4-dimensional 
parent has the form F(T, c, f2) = F(T, c)/47r. We could, in fact, have worked our whole 
procedure in 4 dimensions instead of 2. Had we done so, isotropy of the parent would have 
been imposed at each iteration by explicitly smoothing away all of the VL dependence from 
the stacked Pj(T, c, f2) functions. One would expect, for a plausible set of models leading 
to a plausible parent distribution, that the stacked Pj's should have an VL dependence not 
too far from isotropic, before it is smoothed away. This gives us an important consistency 
check: the sum, SiPj(fi), of the final posterior densities from equation (^) ought to be 
reasonably fiat. Even though the parent distribution is, by construction, isotropic, there is 
no guarantee that the sample is isotropic. If we find a strong orientation bias in the sample 
despite assuming an isotropic parent, this constitutes a contradiction and signals a false 
assumption. 

Two applications of this test are shown in the bottom four panels of Figure Figure 
He shows the parent distribution derived under the assumption that all galaxies rotate 
about their intrinsic long axes (C = 0). Most objects are close to oblate and quite flat, with 
a small but significant fraction of rounder, triaxial systems. This is clearly different from 
the maximal-ignorance distribution in Figure However, this case can be ruled out by 
the orientation distribution of the sample, shown in Figure For the galaxies all to be 
long-axis rotators, we must be seeing them in nearly the same orientation; the line of sight 
lies inside one of two 45 deg-wide cones for about 40% of the sample. 

A better quantitative measure of the orientation bias is the rms deviation of SjPj(fi) 
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from perfect isotropy, normalized to unit mean; we refer to this as the sample anisotropy , Ag. 
For the case in Figure |^e, Ag = 1.17. Table |^ gives Ag values for the distributions calculated 
using various subsets of the dynamical models. Unfortunately, it is not straightforward to 
link a value of Ag with a confidence limit. The expected Ag distribution for an ensemble of 
random isotropic samples depends on the forms of the individual Pj's, which depends on 
both data and models. We can make a very rough correspondence to an easier statistical 
problem if we imagine that each Pi{Q) simply marks a fraction / of the sphere as allowed 
and a fraction 1 — / as excluded. The Pj's are then n patches thrown down at random onto 
the sphere. At a random point on the sphere, the number m of overlapping patches is given 
by a binomial distribution. In Table ^, we find that, except for the C = and C = oo 
cases, all of the models using the kinematic data hover around Ag ^ 0.2. This would follow 
from the binomial distribution for / = 0.63, which is a not-unreasonable characterization 
of the Pj's. We calculate that, if 0.2 is the expected Ag for a random sample of 13 objects 
and if Ag is distributed as in the patch problem, we can reject cases with Ag > 0.40 at 99% 
confidence and cases with Ag > 0.48 at 99.9% confidence. A more realistic simulation with 
9 patches that exclude 50% of the sphere and 4 that exclude only 10% gives very similar 
results. Thus the hypothesis that elliptical galaxies rotate about their long axis is firmly 
ruled out. Of course, this is neither a particularly surprising nor new result; Binney (1985) 
reached the same conclusion from essentially the same data. 

Figure |(i shows the parent distribution under the assumption that all objects rotate 
around their intrinsic short axes (C = oo, also known as "zero intrinsic misalignment"). 
Here, most objects are triaxial, again very different from the maximal ignorance result. 
Figure |^/ shows the orientation distribution for the sample, which has Ag = 0.43. The 
assumption of zero intrinsic misalignment for all systems is excluded at approximately the 
99.6% confidence level. This result differs from that of Franx et al. (1991), who were able 
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to reproduce the observed distribution of ellipticities and kinematic misalignment angles^ 
with a family of triaxial models rotating about their short axes. We have not explored the 
source of this disagreement in depth. While it may be due simply to our smaller sample, 
we suspect that the models with which Franx et al. can fit galaxies with large kinematic 
misalignments fail on more detailed comparison with multi-position-angle data. 



4.3. Dynamical Configurations That Cannot Be Ruled Out 

Of the parent distributions we have derived from various subsets of the dynamical 
models, we find no other cases that can be ruled out on the basis of the Ag values. Some 
of the unexcludable cases nonetheless differ significantly from the maximal ignorance 
distribution. Of particular interest are the cases in the last four rows of Table ^ for which 
the derived parent distributions are shown in Figure ^ These distributions differ only 
in whether the mean rotation is assumed to be disklike or spheroidlike (see § |2.4|) . The 
parent distribution is more sensitive to this assumption than to any of the other dynamical 
parameters, save for the cases already ruled out above. 

Figure ^ shows the triaxiality distributions for these four cases, indicating the 
fraction of nearly oblate, nearly prolate, and triaxial systems as defined in § [4.1| . The 
prevalence of axisymmetric systems over triaxial ones is significantly affected by the rotation 
characteristics, in the sense that axisymmetry becomes less common if rotation is more 
spheroidlike. Moreover, the fractions of nearly prolate and nearly oblate objects are largely 
determined, respectively, by the character of the rotation in the long-axis and short-axis 
tubes. The peak at the prolate limit seen in the maximal ignorance result disappears 

^For galaxies with only major and minor axis kinematics, the misalignment angle is 

tan (^minor/^major) • 
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entirely if the rotation in the long-axis tubes is spheroidhke. The dominant peak at the 
oblate limit is lowered by nearly a factor of two if the short-axis tube rotation is spheroidlike 
rather than disklike. 

We consider this the most important result in this paper: if rotation in ellipticals is 
generally disklike, then triaxiality is rare; if spheroidlike, triaxiality is common. It follows 
that understanding the shapes of elliptical galaxies is closely linked with understanding 
whether weak disks are common structural components. It also follows that a physical 
understanding of what conditions during formation are likely to impose disklike or 
spheroidlike rotation on a hot stellar system would be extremely valuable. 

5. Discussion 

5.1. Previous Results on the Shape Distribution 

A number of attempts have been made in the past to determine the parent intrinsic 
shape distribution of elliptical galaxies, mostly using photometry alone. It is interesting to 
see how our maximal-ignorance result compares to some of these. 

Ryden (1992) fits a parent distribution to a sample of 171 measured ellipticities by 
letting the distribution assume the form of a circular Gaussian in axis ratio space. She 
finds a best-fit center to the distribution at 6 = 0.98, c — 0.69, implying that the most 
common shape is nearly oblate. The distribution is wide, however, with 61% of galaxies 
having a triaxiality between 0.2 and 0.8, compared to 35% for our sample. This difference 
may be attributable to Ryden's assumption of a single peak; our method shows that the 
distribution may be bimodal. When recast in terms of (T, c) , Ryden's distribution has a 
short-to-long axis ratio expectation value (c) = 0.68, similar to our value of 0.71. 

Lambas et al. (1992) take a similar approach using 2135 measurements of ellipticities 
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from the APM Bright Galaxy Survey. Using a Monte Carlo technique they find the 
elliptical Gaussian in axis ratio space which best reproduces their observations. Their 
results are remarkably difi^erent both from ours and from Ryden's. They find the center of 
their distribution at a fiattening c = 0.55 with a width of 0.2 in that dimension, implying 
that 30% of eUipticals have c < 0.4. The main reason for this difference is an excess of 
flat galaxies in their sample. Only 2% of the galaxies in Ryden's (1992) sample have an 
apparent ellipticity e > 0.6, but the APM sample has 30% to 40% in that range. Lambas 
et al. do not offer an explanation for this apparent inconsistency with previous photometric 
studies of eUiptical galaxies. Conceivably a large SO contamination could be the cause. 

A nonparametric, maximum-entropy shape distribution for the Ryden (1992) ellipticity 
sample is derived using a modified Lucy's method by Statler (1994a). He finds a rather 
broad distribution in triaxiality, with 47% of the galaxies having T < 0.5 compared with 
70% in our distribution. 

Using the same data as Statler (1994a), Tremblay and Merritt (1995) use a 
nonparametric maximum penalized likelihood estimator to derive the maximum-entropy 
shape distribution. They find that it is weakly bimodal and weighted towards oblate figures. 
Our distribution is significantly more bimodal and predicts fewer triaxial galaxies. 

Using the same technique, Tremblay and Merritt (1996) estimate the parent distribution 
from a sample of 220 ellipticities. They assume that all galaxies have the same triaxiality 
and then proceed to calculate the distribution of intrinsic fiattenings c. They find that a 
pure oblate or prolate distribution is inconsistent with the available data and that a division 
of intrinsic fiattenings exists between bright and faint galaxies with peaks at c — 0.75 and 
c = 0.65 respectively. All our galaxies are bright and therefore our expectation value of 
c = 0.71 agrees well with theirs. 

Although the above studies, with the exception of Lambas et al. (1992), give similar 
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results for the axis ratio c/a, none is able to put any real constraints on triaxiality, even 
when large samples are used. This demonstrates the need to include kinematic data in the 
models. Franx, lUingworth and de Zeeuw (1991) attempt to address this need by including 
the misalignment between the photometric and kinematic axes in their models. Studying 
a sample of 38 ellipticals, they conclude that a wide variety of distributions are consistent 
with the data, including ones similar to ours with both an oblate and a prolate peak. 



5.2. Previous Results for Individual Galaxies 

Some of the individual galaxies in our sample have been modeled previously. Statler 
(1994c) treats NGC 3379 using essentially the same data and methods applied here, except 
that the galaxy is fit in isolation, using a flat parent distribution. The result is that flattened 
nearly oblate shapes or rounder triaxial configurations are allowed by the data. Compared 
with this earlier result, the posterior density shown in Figure ^jis more constrained toward 
small T due to the preference for near axisymmetry in the parent distribution. 

Some objects in our sample have available additional kinematic or morphological 
constraints which are not included in our models. The best-studied example is NGC 1052, 
which has been modeled by Binney et al. (1990), Tenjes et al. (1993), and Plana and 
Boulesteix (1996). This galaxy has the best constrained shape in our sample; Figure ^ 
shows only a small permitted region around the oblate spheroid with c = 0.63. The small 
triaxiality supports the use of axisymmetric models by Binney et al. (1990) to constrain 
the phase space distribution function. Applying the Jeans equation to the observed surface 
photometry and comparing the predicted velocity dispersion and azimuthal streaming to the 
observed kinematics, they find that NGC 1052 is consistent with a two integral distribution 
function. Tenjes et al. (1993), using the method of Franx, Illingworth and de Zeeuw (1991) 
and the presence of a gas disk to constrain the viewing angles, find that, depending on 
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the specific kinematic model used, c fies between 0.4 and 0.6 and the triaxiality is well 
constrained between 0.56 and 0.61. These values imply a very highly triaxial galaxy, and 
lie well outside of our 95% highest posterior density region. Using similar methods Plana 
& Boulesteix (1996) calculate the triaxiality to be 0.48 with a fiattening of 0.5. This is 
much flatter and more triaxial than our result. It is possible that including orientation 
constraints from the gas would alter our derived shape. However, the results of Tenjes et 
al. (1993) and Plana and Boulesteix (1996) are very sensitive to the orientation of the disk, 
and consequently to assumptions about its intrinsic flatness and circularity; even a small 
error here could change their results dramatically. 

In a study similar to that of Binney et al. (1990), Van der Marel et al. (1990) model 
the distribution functions of NGC 3379, NGC 4261, and NGC 4472. Their use of oblate 
axisymmetric models for NGC 3379 and NGC 4472 is supported by our results for these 
galaxies. For NGC 4261 they flt the observations to a prolate model with c = 0.59, which 
is consistent with our triaxiality estimate and is within our 95% highest posterior density 
region. 

At the risk of disappointing the reader, we have avoided discussing the orientations of 
individual galaxies in the sample. This is, admittedly, counter to the original motivation 
of DB, which was to determine if there is any relationship between the orientations of the 
galaxies and their radio jets. Although we have calculated orientation constraints for each 
of the galaxies, a full discussion of this topic would of necessity be lengthy, and is outside 
the scope of this paper. We will deal with this issue in a future publication. 
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6. Summary and Conclusions 

By combining photometric and kinematic data with dynamical models using the 
method of Statler (1994b), we have derived constraints on the intrinsic shapes and 
orientations of 13 ellipticals from the Davies and Birkinshaw (1988) sample of radio galaxies. 
Using an iterative Bayesian approach we have then combined those results to estimate the 
parent shape distribution from which they were drawn, under the assumption that this 
parent distribution has no preferred orientation. In the process we have obtained improved 
constraints on the shapes of the individual objects. 

We have found that the parent shape distribution shows a tendency toward bimodality, 
with peaks at the oblate and prolate limits. In the distribution derived under minimal 
assumptions about the galaxies' internal dynamics, only about a one-third of the objects 
would be strongly triaxial (0.2 < T < 0.8). However, the parent distribution does depend 
on dynamical assumptions. Some of these assumptions can be ruled out because they would 
require the sample to have a strong orientation bias; configurations in which all galaxies 
rotate purely about either their long axes or their short axes can be excluded on these 
grounds. On the other hand, configurations in which the mean motions in the short-axis and 
long-axis tube orbits are either disklike — dropping off away from the symmetry planes — or 
spheroidlike — staying approximately constant at a given radius — cannot be distinguished at 
this point. Whether the rotation is diskhke or spheroidhke has a strong effect on the inferred 
shape distribution. Spheroidlike rotation in the long-axis or short-axis tubes, respectively, 
significantly reduces the fraction of nearly prolate or nearly oblate galaxies; bimodality 
is completely eliminated if the long-axis tubes are spheroidlike and the short-axis tubes 
disklike. In a nutshell, if rotation in ellipticals is generally disklike, then triaxiality is rare; 
if spheroidlike, triaxiality is common. 

This inferential link between diskiness and axisymmetry complements the intuitive 



-28- 



physical notion that the two ought to go hand in hand. There is evidence from the width 
of the TuUy-Fisher relation that the disks of spiral galaxies are very nearly circular (Pranx 
& de Zeeuw 1992), and indications from numerical experiments that growing even a weak 
disk in a triaxial halo can render the latter axisymmetric (Dubinski 1994). Whether weak 
disks in elliptical galaxies are detectable is another long-standing issue receiving renewed 
attention (Magorrian 1999). High- accuracy, multi-position-angle kinematic mapping may 
be able to reveal hidden disks, but the expected signatures are subtle. Some support 
is lent to the possibility that weak disks may be common by the kinematic similarities 
that the "standard eUiptical" NGC 3379 shares with the SO galaxy NGC 3115 (Statler 
& Smecker-Hane 1999). Theoretically, however, the origin of these particular kinematic 
features is not understood. As we have stressed, a physical understanding of the processes 
that may establish disklike or spheroidlike rotation in a hot stellar system is sorely needed. 

We are indebted to Barbara Ryden and the referee, David Merritt, for numerous 
constructive comments. This work was supported by NASA Astrophysical Theory Program 
Grant NAG5-3050 and NSF CAREER grant AST-9703036. 
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Table 1. Observational data used in models 
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Table 1 — Continued 
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NGC 4636 
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^Degrees from minor axis 



''Only one value of k needed since the logarithmic slope is almost constant between Rmm and Rman 
'^Somces for Kinematic Data: (1) Davies & Birkinshaw (1988) (2) Franx et al. (1989) (3) Binney et al. 
(1990) (4) Bender et al. (1994) (5) Fried & lUingworth (1994) 
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Table 2. Sample Anisotropies from Various Dynamical Assumptions 



Case 


Restriction 




NK 


No kinematics, photometry only 


0.10 





Maximal Ignorance (all models) 


0.18 


1 


C(T,c) = "Prescription 1'"^ 


0.19 


2 


C(T,c) = "Prescription 2'"^ 


0.16 


3 


C(T, c) = "Prescription 3'"^ 


0.20 


4 


C(T, c) "Prescription 4'"^ 


0.17 


5 


C — CO (around short axis) 


0.43 


6 


C = 1 


0.17 


7 


C = 0.5 


0.26 


8 


C = (around long axis) 


1.17 


XI 


sphcroidlikc / sphcroidlike'' 


0.32 


X2 


disklike/ disklike'' 


0.20 


X3 


spheroidlike /disklike'' 


0.21 


X4 


(lisklik(\/spli('r()i<llikc'' 


0.13 



'^Prescriptions 1-4 given by equations (11) - 
(14) of Statler (1994b). 

^Form of v*{t) for long axis/short axis tubes. 
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Fig. 1. — Stack-smooth-iterate algorithm in 1 dimension, finding the distribution of the 
intrinsic quantity X from measurements of an observable x. Likelihoods Li{X) for 9 
measured x values (a) are multiplied by the initial model parent distribution (6), summed 
(c), and smoothed to produce an improved parent {d). Result after 20 iterations {e, solid 
line) is a decent match to the true parent distribution {dotted line). Estimates of X values 
for individual objects (/) differ from the original likelihoods. 

Fig. 2. — {a-d) Parent intrinsic shape distributions, in the space of triaxiality T of the mass 
distribution and flattening of the light distribution, obtained under different assumptions 
for the internal dynamics of the galaxies. In each panel, round prolate galaxies are at upper 
left, flattened oblate galaxies at lower right; objects in between are triaxial. Contours enclose 
68% and 95% of the total probability, (a) "Maximal ignorance" distribution, an unweighted 
combination of all models; (6) result from photometry only, omitting all kinematic data; 
(c) rotation solely around the short axis; (d) rotation solely around the long axis, (e, /) 
Orientation distribution of the sample, for the shape distributions in (c, d) respectively. 
Centers of the circles correspond to views down the short axis, top and bottom edges to 
views down the long axis, and extreme right-hand edge to views down the intermediate axis. 
These cases can be ruled out because of the strong orientation bias implied for the sample. 

Fig. 3. — (a) Triaxiality distribution in the maximal ignorance case, obtained by integrating 
Figure l^a over flattening c. Dashed lines demarcate nearly prolate (left), triaxial (center), 
and nearly oblate (right) regions. Percentages indicate fractions of total probability in each 
region, (b) Expected distribution of apparent ellipticities in the maximal ignorance case 
{smooth curve). Histogram shows ellipticity distribution for a sample of 165 galaxies from 
Ryden (1992). 

Fig. 4. — Posterior probability densities in the intrinsic shape plane for each of the 9 galaxies 
which show signiflcant rotation, using the maximal ignorance parent distribution. Contours 
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indicate the 68% and 95% highest posterior density regions. 

Fig. 5. — Parent intrinsic shape distributions obtained under different assumptions for the 
"diskhke" or "spheroidhke" character of the rotation in long-axis and short-axis tubes, as 
indicated. 

Fig. 6. — Triaxiahty distributions for the same four cases in in fig. ^, obtained by integrating 
the parent distributions over c. Spheroidhke rotation in long-axis or short-axis tubes 
suppresses the fraction of prolate or oblate objects, respectively. 
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